library("rstudioapi")     
##Set working directory
setwd(dirname(getActiveDocumentContext()$path))
getwd() 
library("readxl") 
require(miceadds)
library(data.table)
library(tidyverse)
library(DescTools)
library(htmlTable)
library(stargazer)
library(estimatr)

# Load final Data file: 
load("FinalFinal20062023Full1-6.RData")
load("FinalFinalSP20062023Full1-6.RData")

library(dplyr)
final_finalV2 = mutate(final_finalV2, 
                       Agec = (Age - 30)/10, # per 10 year age
                       Incomec = (Income - 100)/50, # per 50 income
                       clinic_distancec = (clinic_distance - 5)/5, # per 5
                       SNMetricc = SNMetric/10)

final_finalV2P = final_finalV2[which(final_finalV2$individual_treatment == "Placebo"),]

ModP1  <- glm.cluster(vaccine_intention ~ Vcash_cat  + Gender + Agec + Education + Incomec + Employed + SNMetricc + clinic_distancec + 
                        dist2 + dist3 + dist4 + dist5 + dist6, 
                      data = final_finalV2P, family=binomial,
                      cluster = final_finalV2P$Q123)
exp(cbind("Odds ratio" = coef(ModP1), confint.default(ModP1, level = 0.95)))


ModP2 <- glm.cluster(vaccine_reported_combo ~ Vcash_cat  + Gender + Agec + Education + Incomec + Employed + SNMetricc + clinic_distancec + 
                    dist2 + dist3 + dist4 + dist5 + dist6, 
                    data = final_finalV2P, family=binomial,
                    cluster = final_finalV2P$Q123)
exp(cbind("Odds ratio" = coef(ModP2), confint.default(ModP2, level = 0.95)))

ModP3 <- glm.cluster(ActVacApril ~ Vcash_cat  + Gender + Agec + Education + Incomec + Employed + SNMetricc + clinic_distancec + 
                       dist2 + dist3 + dist4 + dist5 + dist6, 
                     data = final_finalV2P, family=binomial,
                     cluster = final_finalV2P$Q123)
exp(cbind("Odds ratio" = coef(ModP3), confint.default(ModP3, level = 0.95)))


# P1 and P2: 
tab_P1 = exp(cbind("Odds ratio" = coef(ModP1), confint.default(ModP1, level = 0.95)))
tab_P2 = exp(cbind("Odds ratio" = coef(ModP2), confint.default(ModP2, level = 0.95)))
tab_P3 = exp(cbind("Odds ratio" = coef(ModP3), confint.default(ModP3, level = 0.95)))


rownames(tab_P1) = c("Intercept", "CashPlacebo", "HealthPlacebo", "Male", "Age (+10 yrs)", "High-Educate", "Low-Educate", "Medium-Educate", "Mean-Food (+50)", "Employed", "Social Media (+10)", "Access (+5)","District 2", "District 3", "District 4", "District 5","District 6")
rownames(tab_P2) = c("Intercept", "CashPlacebo", "HealthPlacebo", "Male", "Age (+10 yrs)", "High-Educate", "Low-Educate", "Medium-Educate", "Mean-Food (+50)", "Employed", "Social Media (+10)", "Access (+5)","District 2", "District 3", "District 4", "District 5","District 6")
rownames(tab_P3) = c("Intercept", "CashPlacebo", "HealthPlacebo", "Male", "Age (+10 yrs)", "High-Educate", "Low-Educate", "Medium-Educate", "Mean-Food (+50)", "Employed", "Social Media (+10)", "Access (+5)","District 2", "District 3", "District 4", "District 5","District 6")


convert_fun = function(data, name_mod){
  new_vec = c()
  for (i in c(1:dim(data)[1])){
    new_vec[i] = as.character(paste(format(round(data[i,1],2), nsmall = 2),paste("(",paste(format(round(data[i,2],2), nsmall = 2), format(round(data[i,3],2), nsmall = 2), sep = ", "),")", sep ="")))
  }
  new_vec = data.frame(new_vec)
  rownames(new_vec) = rownames(data)
  colnames(new_vec) = name_mod
  return(new_vec)
}

col_P1 = convert_fun(tab_P1, "Vaccine Intention")
col_P2 = convert_fun(tab_P2, "Vaccine Reported")
col_P3 = convert_fun(tab_P2, "Actual Vaccination")

tb1 = merge(col_P1, col_P2, by = 0, all = TRUE);rownames(tb1) = tb1[,1];tb1[,1] = NULL
tb2 = merge(tb1, col_P3, by = 0, all = TRUE);rownames(tb2) = tb2[,1];tb2[,1] = NULL

no.obs = c((ModP1$glm_res$df.null),(ModP2$glm_res$df.null), (ModP3$glm_res$df.null))
aic_mod = round(c(ModP1$glm_res$aic, ModP2$glm_res$aic, ModP3$glm_res$aic))
log_lik = round(c(as.vector(logLik(ModP1$glm_res)), as.vector(logLik(ModP2$glm_res)), as.vector(logLik(ModP3$glm_res))))
mod_spec = data.frame(rbind(no.obs, aic_mod, log_lik))
rownames(mod_spec) = c("Observations", "Akaike Inf. Crit.", "Log Likelihood")
colnames(mod_spec) = colnames(tb2)

tab_P = rbind(tb2, mod_spec)
tab_P2 = tab_P[c("CashPlacebo","HealthPlacebo", "Access (+5)" ,"Age (+10 yrs)", "Male", "High-Educate", "Medium-Educate", "Low-Educate","Employed","Mean-Food (+50)","Social Media (+10)","District 2", "District 3", "District 4", "District 5","District 6", "Intercept", "Observations", "Akaike Inf. Crit.", "Log Likelihood"),]
stargazer(tab_P2, type = "text", summary = FALSE, out="SM_TablesS2.tex")
stargazer(tab_P2, type = "latex", summary = FALSE, out="SM_TablesS2.tex")
stargazer(tab_P2, type = "latex", summary = FALSE)








